// project home: cp3-cg.googlecode.com

#include "myvectormath.h"
#include "global.h"
#include <iostream>
#include <cmath>
using namespace std;

vector<Double_t> matrixvectorproduct(vector<Double_t> *A, vector<Int_t> *JA, vector<Int_t> *IA, vector<Double_t> *x, bool debug){
  
//   bool debug=true;
  
  Int_t dim = x->size();
  if (debug){cout << "dimension: " << dim << endl;}
  vector<Double_t> b(dim,0); //Ausgabevektor
  
  for (int j=0;j<dim;j++){
    if (debug){cout << "j=" << j << ";" << endl;}
    for (int i=IA->at(j);i<IA->at(j+1);i++){
      if (debug){cout << "i=" << i << ";" << endl;}
      b.at(j)=b.at(j)+A->at(i)*x->at(JA->at(i));
    }
  }
  
  return b;
}

void Laplaceproduct(Int_t n, vector<Double_t>  *vecin, vector<Double_t>  *vecout){
  /*
   * Angelegt und berechnet wird
   * Volumen       nvol
   * NN-Indexfeld  nn[k][i], k=0..2*ndim, i=0..(nvol-1)
   * 
   * nn[k][i] gibt den Index des Nachbarn von i in Richtung +k,
   * unter Beruecksichtigung periodischer Randbedingungen.
   * Fuer einen Schritt in Richtung -k setze man den Index auf (ndim+k).
   * nn[0][i] ist reserviert.*/
  
  
  // nn[][] nachbarn
  // ndim dimension
  
  for (int i=0;(unsigned)i<vecout->size();i++){
    vecout->at(i)=(2*ndim+pow(mass,2))*vecin->at(i);
    for (int k=1;k<=ndim;k++){
      vecout->at(i)=vecout->at(i)-(vecin->at(nn[k][i])+vecin->at(nn[k+ndim][i]));
    }
  }
  
}








Double_t skalarp(vector<Double_t> *x, vector<Double_t> *y){
  Double_t product=0;
  for (int i=0; (unsigned)i<x->size();i++){
    product=product+x->at(i) * y->at(i);
  }
  return product;
}


void vectorsub(vector<Double_t> *x, vector<Double_t> *y, vector<Double_t> *vecout){
  vecout->assign(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    vecout->at(i)=x->at(i) - y->at(i);
  }
}


void vectoradd(vector<Double_t> *x, vector<Double_t> *y, vector<Double_t> *vecout){
  vecout->assign(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    vecout->at(i)=x->at(i) + y->at(i);
  }
}


void VectorMultiS(vector<Double_t> *x, Double_t alpha, vector<Double_t> *vecout){
  vecout->assign(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    vecout->at(i)=x->at(i) * alpha;
  }
}





